function [im,snrP,imS] = textureSynthesis15D(params, im0, Niter, cmask, imask)

% [res,snrP,imS] = textureSynthesis15D(params, initialIm, Niter, cmask, imask)
%
% Synthesize 1.5D textures by applying Portilla-Simoncelli model/algorithm.
%
% params: structure containing texture parameters (as returned by textureAnalysis15D).
%
% im0: initial image, OR a vector (Ydim, Xdim, [SEED]) containing
% dimensions of desired image and an optional seed for the random
% number generator.  If dimensions are passed, initial image is
% Gaussian white noise.
%
% Niter (optional): Number of iterations.  Default = 50.
%
% cmask (optional): binary column vector (4x1) indicating which sets of
% constraints we want to apply in the synthesis. The four sets are:
%               1) Marginal statistics (mean, var, skew, kurt, range)
%               2) Correlation of subbands (space, orientation, scale)
%               3) Correlation of magnitude responses (sp, or, sc)
%               4) Relative local phase
%
% imask (optional): imsizex2 matrix.  First column is a mask, second
% 	column contains the image values to be imposed. If only one column is
%	provided, it assumes it corresponds to the image values, and it uses
%	a raised cosine square for the mask.
% snrP (optional):	Set of adjustment values (in dB) of the parameters.
% imS (optional):	Sequence of synthetic images, from niter = 1 to 2^n, being
% 			n = floor(log2(Niter)).

% Modifications for 1.5D textures are described in:
%  "Self-Similar Texture for Coherent Line Stylization,"
%  Pierre Benard, Forrester Cole, Aleksey Golovinskiy, Adam Finkelstein,
%  NPAR 2010.
%
% Original code by:
% Javier Portilla and Eero Simoncelli.
% Work described in:
%  "A Parametric Texture Model based on Joint Statistics of Complex Wavelet Coefficients".
%  J Portilla and E P Simoncelli. Int'l Journal of Computer Vision,
%  vol.40(1), pp. 49-71, Dec 2000.
%
% Please refer to the original publication if you use the program for research or
% for technical applications. Thank you.
%
% Copyright, Center for Neural Science, New York University, January 2001.
% All rights reserved.

Warn = 0;  % Set to 1 if you want to see warning messages

%% Check required args are passed:
if (nargin < 2)
    error('Function called with too few input arguments');
end

if ( ~exist('Niter') | isempty(Niter) )
    Niter = 50;
end

if (exist('cmask') & ~isempty(cmask) )
    cmask = (cmask > 0.5);  % indices of ones in mask
else
    cmask = ones(4,1);
end

%% Extract parameters
statg0 = params.pixelStats;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Grab the horizontal statistics.
% -fcole May 17 2010

horMeans0 = params.horMeans0;
horVars0 = params.horVars0;

% End modification
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mean0 = statg0(1); var0 =  statg0(2);
skew0 =  statg0(3); kurt0 =  statg0(4);
mn0 =  statg0(5);  mx0 = statg0(6);
statsLPim = params.pixelLPStats;
skew0p = statsLPim(:,1);
kurt0p = statsLPim(:,2);
vHPR0 = params.varianceHPR;
acr0 = params.autoCorrReal;
ace0 = params.autoCorrMag;
magMeans0 = params.magMeans;
C0 = params.cousinMagCorr;
Cx0 = params.parentMagCorr;
Crx0 = params.parentRealCorr;

%% Extract {Nsc, Nor, Na} from params
tmp = size(params.autoCorrMag);
Na = tmp(1); Nsc = tmp(3);
Nor = tmp(length(tmp))*(length(tmp)==4) + (length(tmp)<4);
la = (Na-1)/2;

%% If im0 is a vector of length 2, create Gaussian white noise image of this
%% size, with desired pixel mean and variance.  If vector length is
%% 3,  use the 3rd element to seed the random number generator.
if ( length(im0) <= 3 )
    if ( length(im0) == 3)
        randn('state', im0(3)); % Reset Seed
        im0 = im0(1:2);
    end
    im = mean0 + sqrt(var0)*randn(im0);
else
    im = im0;
end

%% If the spatial neighborhood Na is too big for the lower scales,
%% "modacor22.m" will make it as big as the spatial support at
%% each scale:
[Ny,Nx] = size(im);
nth = log2(min(Ny,Nx)/Na);
if nth<Nsc+1 & Warn,
    fprintf(1,'Warning: Na will be cut off for levels above #%d !\n',floor(nth));
end

if  exist('imask') & ~isempty(imask),
    if  size(imask,1) ~= prod(size(im))
        error(sprintf('imask size %d does not match image dimensions [%d,%d]',...
            size(imask,1), size(im,1), size(im,2)));
    end
    if size(imask,2) == 1,
        mask = (cos(-pi/2:2*pi/Ny:pi*(1-2/Ny)/2)).'*cos(-pi/2:2*pi/Nx:pi*(1-2/Nx)/2);
        mask = mask.^2;
        aux = zeros(size(im));
        aux(Ny/4+1:Ny/4+Ny/2,Nx/4+1:Nx/4+Nx/2)=mask;
        mask = aux;
    else
        mask = reshape(imask(:,1),size(im));
    end
end

%imf = max(1,gcf-1); snrf = imf+1;
imf = 555; snrf = imf+1;
figure(imf);  clf
subplot(1,2,1); grayRange = showIm(im,'auto',1); title('Starting image');
drawnow

prev_im=im;

%imf = max(1,gcf-1);
imf = 555;
figure(imf);
clf;showIm(im,'auto',1); title(sprintf('iteration 0'));

nq = 0;
Nq = floor(log2(Niter));
imS = zeros(Ny,Nx,Nq);

%% MAIN LOOP
for niter = 1:Niter

    %p = niter/Niter;
    p = 1;

    %% Build the steerable pyramid
    [pyr,pind] = buildSCFpyr(im,Nsc,Nor-1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Commented out the following error.
% -fcole May 17 2010

    %if ( any(vector(mod(pind,2))) )
        %error('Algorithm will fail: band dimensions are not all multiples of 4!');
    %end
% End modification
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %% Subtract mean of lowBand:
    nband = size(pind,1);
    pyr(pyrBandIndices(pind,nband)) = ...
        pyrBand(pyr,pind,nband) - mean2(pyrBand(pyr,pind,nband));

    apyr = abs(pyr);

    %% Adjust autoCorr of lowBand
    nband = size(pind,1);
    ch = pyrBand(pyr,pind,nband);
    Sch = min(size(ch)/2);
    nz = sum(sum(~isnan(acr0(:,:,Nsc+1))));
    lz = (sqrt(nz)-1)/2;
    le = min(Sch/2-1,lz);
    im = real(ch);  %Reconstructed image: initialize to lowband
    [mpyr,mpind] = buildSFpyr(im,0,0);
    im = pyrBand(mpyr,mpind,2);
    vari =  acr0(la+1:la+1,la+1:la+1,Nsc+1);

    if cmask(2),
        if vari/var0 > 1e-4,
            [im, snr2(niter,Nsc+1)] = ...
                modacor22(im, acr0(la-le+1:la+le+1,la-le+1:la+le+1,Nsc+1),p);
        else
            im = im*sqrt(vari/var2(im));
        end
        if (var2(imag(ch))/var2(real(ch)) > 1e-6)
            fprintf(1,'Discarding non-trivial imaginary part, lowPass autoCorr!');
        end
        im = real(im);
    end % cmask(2)
    if cmask(1),
        if vari/var0 > 1e-4,
            [im,snr7(niter,2*(Nsc+1)-1)] = modskew(im,skew0p(Nsc+1),p);	% Adjusts skewness
            [im,snr7(niter,2*(Nsc+1))] = modkurt(im,kurt0p(Nsc+1),p);	% Adjusts kurtosis
        end
    end	% cmask(2)

    %% Subtract mean of magnitude
    if cmask(3),
        magMeans = zeros(size(pind,1), 1);
        for nband = 1:size(pind,1)
            indices = pyrBandIndices(pind,nband);
            magMeans(nband) = mean2(apyr(indices));
            apyr(indices) = apyr(indices) - magMeans(nband);
        end
    end	% cmask(3)

    %% Coarse-to-fine loop:
    for nsc = Nsc:-1:1

        firstBnum = (nsc-1)*Nor+2;
        cousinSz = prod(pind(firstBnum,:));
        ind = pyrBandIndices(pind,firstBnum);
        cousinInd = ind(1) + [0:Nor*cousinSz-1];

        %% Interpolate parents
        if 1%,(cmask(3) | cmask(4)),
            if (nsc<Nsc)
                parents = zeros(cousinSz,Nor);
                rparents = zeros(cousinSz,Nor*2);
                for nor = 1:Nor
                    nband = (nsc+1-1)*Nor+nor+1;

                    tmp = expand(pyrBand(pyr, pind, nband),2)/4;
                    rtmp = real(tmp); itmp = imag(tmp);
                    tmp = sqrt(rtmp.^2 + itmp.^2) .* exp(2 * sqrt(-1) * atan2(rtmp,itmp));
                    rparents(:,nor) = vector(real(tmp));
                    rparents(:,Nor+nor) = vector(imag(tmp));

                    tmp = abs(tmp);
                    parents(:,nor) = vector(tmp - mean2(tmp));
                end
            else
                rparents = [];
                parents = [];
            end
        end % if (cmask(3) | cmask(4))

        if cmask(3),
            %% Adjust cross-correlation with MAGNITUDES at other orientations/scales:
            cousins = reshape(apyr(cousinInd), [cousinSz Nor]);
            nc = size(cousins,2);   np = size(parents,2);
            if (np == 0)
                [cousins, snr3(niter,nsc)] = adjustCorr1s(cousins, C0(1:nc,1:nc,nsc), 2, p);
            else
                [cousins, snr3(niter,nsc), snr4(niter,nsc)] = ...
                    adjustCorr2s(cousins, C0(1:nc,1:nc,nsc), parents, Cx0(1:nc,1:np,nsc), 3, p);
            end
            if (var2(imag(cousins))/var2(real(cousins)) > 1e-6)
                fprintf(1,'Non-trivial imaginary part, mag crossCorr, lev=%d!\n',nsc);
            else
                cousins = real(cousins);
                ind = cousinInd;
                apyr(ind) = vector(cousins);
            end

            %% Adjust autoCorr of mag responses
            nband = (nsc-1)*Nor+2;
            Sch = min(pind(nband,:)/2);
            nz = sum(sum(~isnan(ace0(:,:,nsc,1))));
            lz = (sqrt(nz)-1)/2;
            le = min(Sch/2-1,lz);
            for nor = 1:Nor,
                nband = (nsc-1)*Nor+nor+1;
                ch = pyrBand(apyr,pind,nband);
                [ch, snr1(niter,nband-1)] = modacor22(ch,...
                    ace0(la-le+1:la+le+1,la-le+1:la+le+1,nsc,nor), p);
                ch = real(ch);
                ind = pyrBandIndices(pind,nband);
                apyr(ind) = ch;
                %% Impose magnitude:
                mag = apyr(ind) + magMeans0(nband);
                mag = mag .* (mag>0);
                pyr(ind) = pyr(ind) .* (mag./(abs(pyr(ind))+(abs(pyr(ind))<eps)));
            end
        end   % cmask(3)

        %% Adjust cross-correlation of REAL PARTS at other orientations/scales:
        cousins = reshape(real(pyr(cousinInd)), [cousinSz Nor]);
        Nrc = size(cousins,2);  Nrp = size(rparents,2);
        if  cmask(4) & (Nrp ~= 0)
            a3 = 0; a4 = 0;
            for nrc = 1:Nrc,
                cou = cousins(:,nrc);
                [cou, s3, s4] = ...
                    adjustCorr2s(cou,mean(cou.^2),rparents,Crx0(nrc,1:Nrp,nsc), 3);
                a3 = s3 + a3;
                a4 = s4 + a4;
                cousins(:,nrc) = cou;
            end
            snr4r(niter,nsc) = a4/Nrc;
        end
        if (var2(imag(cousins))/var2(real(cousins)) > 1e-6)
            fprintf(1,'Non-trivial imaginary part, real crossCorr, lev=%d!\n',nsc);
        else
            %%% NOTE: THIS SETS REAL PART ONLY - signal is now NONANALYTIC!
            pyr(cousinInd) = vector(cousins(1:Nor*cousinSz));
        end

        %% Re-create analytic subbands
        dims = pind(firstBnum,:);
        ctr = ceil((dims+0.5)/2);
        ang = mkAngle(dims, 0, ctr);
        ang(ctr(1),ctr(2)) = -pi/2;
        for nor = 1:Nor,
            nband = (nsc-1)*Nor+nor+1;
            ind = pyrBandIndices(pind,nband);
            ch = pyrBand(pyr, pind, nband);
            ang0 = pi*(nor-1)/Nor;
            xang = mod(ang-ang0+pi, 2*pi) - pi;
            amask = 2*(abs(xang) < pi/2) + (abs(xang) == pi/2);
            amask(ctr(1),ctr(2)) = 1;
            amask(:,1) = 1;
            amask(1,:) = 1;
            amask = fftshift(amask);
            ch = ifft2(amask.*fft2(ch));	% "Analytic" version
            pyr(ind) = ch;
        end

        %% Combine ori bands
        bandNums = [1:Nor] + (nsc-1)*Nor+1;  %ori bands only
        ind1 = pyrBandIndices(pind, bandNums(1));
        indN = pyrBandIndices(pind, bandNums(Nor));
        bandInds = [ind1(1):indN(length(indN))];
        %% Make fake pyramid, containing dummy hi, ori, lo
        fakePind = pind([bandNums(1), bandNums, bandNums(Nor)+1],:);
        fakePyr = [zeros(prod(fakePind(1,:)),1);...
            real(pyr(bandInds)); zeros(prod(fakePind(size(fakePind,1),:)),1)];
        ch = reconSFpyr(fakePyr, fakePind, [1]);     % recon ori bands only
        im = real(expand(im,2))/4;
        im = im + ch;
        vari =  acr0(la+1:la+1,la+1:la+1,nsc);
        if cmask(2),
            if vari/var0 > 1e-4,
                [im, snr2(niter,nsc)] = ...
                    modacor22(im, acr0(la-le+1:la+le+1,la-le+1:la+le+1,nsc), p);
            else
                im = im*sqrt(vari/var2(im));
            end
        end	% cmask(2)
        im = real(im);

        if cmask(1),
            %% Fix marginal stats
            if vari/var0 > 1e-4,
                [im,snr7(niter,2*nsc-1)] = modskew(im,skew0p(nsc),p);       % Adjusts skewness
                [im,snr7(niter,2*nsc)] = modkurt(im,kurt0p(nsc),p);         % Adjusts kurtosis
            end
        end	% cmask(1)

    end  %END Coarse-to-fine loop

    %% Adjust variance in HP, if higher than desired
    if (cmask(2)|cmask(3)|cmask(4)),
        ind = pyrBandIndices(pind,1);
        ch = pyr(ind);
        vHPR = mean2(ch.^2);
        if vHPR > vHPR0,
            ch = ch * sqrt(vHPR0/vHPR);
            pyr(ind) = ch;
        end
    end % cmask
    im = im + reconSFpyr(real(pyr), pind, [0]);  %recon hi only

    %% Pixel statistics
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Set the horizontal means and variances 
%   (Section 3.2.2 in "Self-Similar Texture for Coherent Line Stylization")
% -fcole May 17 2010
    
    horMeans = mean(im, 2);
    means = mean2(im);
    horVars = var(im, 1, 2);
    vars = var2(im, means);
    snr7(niter,2*(Nsc+1)+1) = snr(var0,var0-vars);
    for i = 1:size(im, 1)
        im(i, :) = im(i, :) - horMeans(i);  			% Adjusts mean and variance
    end
    [mns mxs] = range2(im + mean0);
    snr7(niter,2*(Nsc+1)+2) = snr(mx0-mn0,sqrt((mx0-mxs)^2+(mn0-mns)^2));
    if cmask(1),
        for i = 1:size(im, 1)
            im(i, :) = im(i, :)*sqrt(((1-p)*horVars(i) + p*horVars0(i))/horVars(i));
        end
    end	% cmaks(1)
    for i = 1:size(im, 1)
        im(i, :) = im(i, :) + horMeans0(i);  			% Adjusts mean and variance
    end

% End modification
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    if cmask(1),
        [im, snr7(niter,2*(Nsc+1)+3)] = modskew(im,skew0,p);	% Adjusts skewness (keep mean and variance)
        [im, snr7(niter,2*(Nsc+1)+4)] = modkurt(im,kurt0,p);	% Adjusts kurtosis (keep mean and variance,
        % but not skewness)
        im = max(min(im,(1-p)*max(max(im))+p*mx0),...
            (1-p)*min(min(im))+p*mn0);		% Adjusts range (affects everything)
    else
        snr7(niter,2*(Nsc+1)+3) = snr(skew0,skew0-skew2(im));
        snr7(niter,2*(Nsc+1)+4) = snr(kurt0,kurt0-kurt2(im));
    end	% cmask(1)

    %% Force pixels specified by image mask
    if (exist('imask') & ~isempty(imask) )
        im = mask.*reshape(imask(:,2 - (size(imask,2)==1)),size(im)) + ...
            (1-mask).*im;
    end

    snr6(niter,1) = snr(im-mean0,im-prev_im);

    if floor(log2(niter))==log2(niter),
        nq = nq + 1;
        imS(:,:,nq) = im;
    end

    tmp = prev_im;
    prev_im=im;

    figure(imf);
    subplot(1,2,1);
    showIm(im-tmp,'auto',1); title('Change');
    subplot(1,2,2);
    showIm(im,'auto',1); title(sprintf('iteration %d/%d',niter,Niter));
    drawnow


    % accelerator
    alpha = 0.8;
    im = im + alpha*(im - tmp);

    commented = 1;  % set it to 0 for displaying convergence of parameters in SNR (dB)
    if ~commented,

        % The graphs that appear reflect
        % the relative distance of each parameter or group
        % of parametersi, to the original's, in decibels.
        % Note, however, that when the original parameters
        % are close to zero, this measurement is meaningless.
        % This is why in some cases it seems that some of
        % the parameters do not converge at all.

        figure(snrf);
        if cmask(1)
            subplot(171); plot(snr7); title('Mrgl stats');
        end
        if cmask(2),
            subplot(172); plot(snr2); title('Raw auto');
        end
        if cmask(3),
            subplot(173); plot(snr1); title('Mag auto');
            subplot(174); plot(snr3); title('Mag ori');
            subplot(175); plot(snr4); title('Mag scale');
        end
        if (Nrp > 0) & cmask(4),
            subplot(176); plot(snr4r); title('Phs scale');
        end
        subplot(177); plot(snr6); title('Im change');
        drawnow

    end  % if ~commented

end %END  MAIN LOOP

im = prev_im;

snrP = [];%snr7 snr2 snr1 snr3 snr4 snr4r snr6];

